#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <mpi.h>
#include <omp.h>

#define NTHREADS 2

int linhas(double **, int);
int colunas(double **, int);

int main(int argc, char *argv[]) {
  
  /*
    nproc - número de processos utilizados
    rank - identificador de cada processo no programa
  */
  int nproc, rank, nescravos;
  MPI_Status status;

  // Inicia execução da MPI
  if(MPI_Init(&argc, &argv) != MPI_SUCCESS) {
    printf("ERRO: Erro ao iniciar MPI.\n");
    MPI_Abort(MPI_COMM_WORLD, -1);
  }

  omp_set_num_threads(NTHREADS);

  MPI_Comm_size(MPI_COMM_WORLD, &nproc); // Recupera o número de processos que o programa está sendo dividido
  MPI_Comm_rank(MPI_COMM_WORLD, &rank); // Recupera o identificador do processo sendo rodado
  nescravos = nproc-1;
	
  // Se for menos de 1 processo (impossível), o programa termina
  if(nproc < 1) {
    printf("ERRO: Numero de processo tem que ser maior que 1.\n");
    MPI_Abort(MPI_COMM_WORLD, -1);
  }
		
  // Se for o processo mestre
  if (rank == 0) {
	
    int J_ORDER, J_ROW_TEST, J_ITE_MAX;
    double J_ERROR, maxDif, maxX, mr=0, sum;
    double **MA, *MB, *x, *xOld, *dif;
    int i, j, escravo, it; // contadores;	
    int *nlinhas, offset, continua=1;
    FILE *f;


    if(argc > 1 && argv[1] != NULL) 
      f = fopen(argv[1], "r");
    else {
      printf("ERRO: Erro ao abrir o arquivo.\n");
      MPI_Abort(MPI_COMM_WORLD, -1);
    }

    // Le o arquivo
    fscanf(f,"%d",&J_ORDER);
    fscanf(f,"%d",&J_ROW_TEST);
    fscanf(f,"%lf",&J_ERROR);
    fscanf(f,"%d",&J_ITE_MAX);

    // Alocações	
    x = (double *) malloc(J_ORDER*sizeof(double));
    xOld = (double *) malloc(J_ORDER*sizeof(double));
    dif = (double *) malloc(J_ORDER*sizeof(double));
	MA = (double **) malloc(J_ORDER*sizeof(double*));
	for(i = 0; i < J_ORDER; i++)
		MA[i] = (double *) malloc(J_ORDER*sizeof(double));
    MB = (double *) malloc(J_ORDER*sizeof(double));
    nlinhas = (int *) malloc(nescravos*sizeof(int));

    for(i = 0; i < J_ORDER; i++)
      for(j = 0; j < J_ORDER; j++)
        fscanf(f, "%lf", &MA[i][j]);

    for(i = 0; i < J_ORDER; i++)
      fscanf(f, "%lf", &MB[i]);

    fclose(f);

    // Verifica a convergência da Matriz A
    if(!linhas(MA, J_ORDER) && !colunas(MA, J_ORDER)) {
      printf("ERRO: Matriz nao converge, portanto o metodo nao funciona.\n");
      MPI_Abort(MPI_COMM_WORLD, -1);
    }

    // Procura o número de linhas que serão dividas nos nós
    offset = ceil((double)J_ORDER/nescravos);

    // Inicia X
    #pragma omp parallel for    
    for(i = 0; i < J_ORDER; i++)
      x[i] = MB[i]/MA[i][i];

    for(escravo=1 ; escravo <= nescravos ; escravo++) {
      if(escravo != nescravos)
        nlinhas[escravo] = offset;
      else
        nlinhas[escravo] = J_ORDER - offset*(nescravos-1); // Se for o ultimo nó, pega o resto das linhas

      // Envia a ordem, o offset e quantas linhas aquela matriz tem
      MPI_Send(&J_ORDER, 1, MPI_INT, escravo, 1, MPI_COMM_WORLD);
      MPI_Send(&offset, 1, MPI_INT, escravo, 1, MPI_COMM_WORLD);
      MPI_Send(&nlinhas[escravo], 1, MPI_INT, escravo, 1, MPI_COMM_WORLD);

      // Envia parte da matriz A, e a matriz B
      for(i = (escravo-1)*offset ; i < (escravo-1)*offset + nlinhas[escravo] ; i++)
        MPI_Send(MA[i], J_ORDER, MPI_DOUBLE, escravo, 1, MPI_COMM_WORLD);
      MPI_Send(MB, J_ORDER, MPI_DOUBLE, escravo, 1, MPI_COMM_WORLD);
    }

    // Trata das iterações e do erro
    for(it=0 ; continua ; it++) {
      for(i = 0; i < J_ORDER; i++)
        xOld[i] = x[i];
      
      // Envia para todos os nós escravos o vetor X da iteração anterior
      MPI_Bcast(xOld, J_ORDER, MPI_DOUBLE, 0, MPI_COMM_WORLD);

      // Recebe os pedaços do novo X calculado nos nós escravos
      for(escravo=1 ; escravo <= nescravos ; escravo++)
        MPI_Recv(&x[(escravo-1)*offset], nlinhas[escravo], MPI_DOUBLE, escravo, 1, MPI_COMM_WORLD, &status);

      // Calcula o erro entre o novo X e o X da iteração anterior
      for(i = 0, maxDif = 0, maxX = 0; i < J_ORDER; i++) {
        sum = x[i] - xOld[i];
        sum = (sum < 0) ? (-1)*sum : sum;
				
        dif[i] = sum;
        maxDif = (sum > maxDif) ? sum : maxDif;

        sum = x[i];
        sum = (sum < 0) ? (-1)*sum : sum;
        maxX = (sum > maxX) ? sum : maxX;
      }

      mr = (maxX != 0) ? maxDif/maxX : mr;

      // Verifica se as iterações acabaram, e manda esse sinal para todos os nós
      continua = (it < J_ITE_MAX && mr > J_ERROR);
      MPI_Bcast(&continua, 1, MPI_INT, 0, MPI_COMM_WORLD);
    }

    // Faz a verificação do algoritmo, para apresentar na tela.
    sum = 0;
    #pragma omp parallel for reduction(+:sum)
    for(i = 0; i < J_ORDER; i++)
      sum = sum + MA[J_ROW_TEST][i]*x[i];

    // Imprime na tela o numero de iterações utilizadas e o teste para verificar se o algoritmo funcionou.
    printf("----------------------------------------\n");
    printf("Iterations: %d\n", it);
    printf("RowTest: %d => [%lf] =? %lf\n", J_ROW_TEST, sum, MB[J_ROW_TEST]);
    printf("----------------------------------------\n");

  } else {
	
    int offset, nlinhas, i, j, k, J_ORDER, continua=1;
    double **MA, *MB; // MA*x = MB
    double *xOld, *partX, sum;

    // Recebe a ordem, offset e número de linhas que tal nó escravo tratará
    MPI_Recv(&J_ORDER, 1, MPI_INT, 0, 1, MPI_COMM_WORLD, &status);
    MPI_Recv(&offset, 1, MPI_INT, 0, 1, MPI_COMM_WORLD, &status);
    MPI_Recv(&nlinhas, 1, MPI_INT, 0, 1, MPI_COMM_WORLD, &status);

    // Aloca espaco para matrizes e vetores.
	MA = (double **) malloc(J_ORDER*sizeof(double*));
	for(i = 0; i < nlinhas; i++)
		MA[i] = (double *) malloc(J_ORDER*sizeof(double));
    MB = (double *) malloc(J_ORDER*sizeof(double));
    xOld = (double *) malloc(J_ORDER*sizeof(double));
    partX = (double *) malloc(nlinhas*sizeof(double));

    // Recebe o pedaço da matriz que esse nó tratará, e recebe o vetor B
    for(i = 0; i < nlinhas; i++)
      MPI_Recv(MA[i], J_ORDER, MPI_DOUBLE, 0, 1, MPI_COMM_WORLD, &status);
    MPI_Recv(MB, J_ORDER, MPI_DOUBLE, 0, 1, MPI_COMM_WORLD, &status);

    // Fica no while, conversando com o nó mestre, até que se ache o melhor resultado
    while(continua) {

      // Recebe do nó pai o X da iteração anterior, para que se possa calcular o novo pedaço de x, partX
      MPI_Bcast(xOld, J_ORDER, MPI_DOUBLE, 0, MPI_COMM_WORLD);

      for(i = (rank-1)*offset, k=0; k < nlinhas; i++, k++) {
        sum = 0;
        //#pragma omp parallel for reduction(+:sum)
        for(j = 0; j < J_ORDER; j++)
          if(i != j)
            sum += MA[k][j]*xOld[j];
        partX[k] = (MB[i] - sum)/MA[k][i];
      }

      // Envia a parte do vetor X calculado, e espera o sinal do mestre para ver se continua ou não no while
      MPI_Send(partX, nlinhas, MPI_DOUBLE, 0, 1, MPI_COMM_WORLD);
      MPI_Bcast(&continua, 1, MPI_INT, 0, MPI_COMM_WORLD);
    }
  }
	
  MPI_Finalize();
  return 0;
}

int linhas(double **MA, int J_ORDER) {
  int i, j;
  double sum, aux, max = 0;

  for(i = 0; i < J_ORDER ; i++) {
    sum = 0;
    if(MA[i][i] == 0)
      return 0;

    for(j = 0; j < J_ORDER; j++)
      if(i != j) {
        aux = MA[i][j]/MA[i][i];
        aux = (aux < 0) ? (-1*aux) : aux;
        sum = sum + aux;
      }

    max = (sum > max) ? sum : max;
  }

  return (max < 1);
}

int colunas(double **MA, int J_ORDER) {
  int i, j;
  double sum, aux, max = 0;

  for(j = 0; j < J_ORDER ; j++) {
    sum = 0;
      if(MA[j][j] == 0)
        return 0;

    for(i = 0; i < J_ORDER; i++)
      if(i != j) {
        aux = MA[i][j]/MA[i][i];
        aux = (aux < 0) ? (-1*aux) : aux;
        sum = sum + aux;
      }

    max = (sum > max) ? sum : max;
  }

  return (max < 1);
}
